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In copper- and iron-based unconventional superconductors, the Bogoliubov quasiparticles interact with a 
spin resonance at momentum (tt, n). This interaction is revealed by specific signatures in the quasiparticle 
spectroscopies, like kinks in photoemission and dips in tunneling. We study these signatures, as they appear 
inside and around a vortex core in the local density of states (LDOS), a property accessible experimentally 
by scanning tunneling spectroscopy. Our model retains the whole nonlocal structure of the self-energy in 
space and time and is therefore not amenable to a Hamiltonian treatment using Bogoliubov-de Gennes 
equations. The interaction with the spin resonance does not suppress the zero-bias peak at the vortex 
center, although it reduces its spectral weight; neither does it smear out the vortex LDOS, but rather it adds 
structure to it. Some of the signatures we find may have been already measured in FeSe, but remained 
unnoticed. We compare the LDOS as a function of both energy and position with and without coupling to 
the spin resonance and observe, in particular, that the quasiparticle interference patterns around the vortex 
are strongly damped by the coupling. We study in detail the transfer of spectral weight induced both locally 
and globally by the interaction and also by the formation of the vortex. Finally, we introduce a new way of 
imaging the quasiparticles in real space, which combines locality and momentum-space sensitivity. This 
approach allows one to access quasiparticle properties that are not contained in the LDOS. 

PACS numbers: 74.55.+V, 74.72.-h, 74.20.Pq 


1. INTRODUCTION 

A renewed interest in topological states of matter^ has 
contributed to bring again into focus the zero modes, which 
appear at the interface between regions with mismatched 
topologies. The electronic states bound to the core of vor¬ 
tices in type-II superconductors belong to this family.^"^ 
They are tied to the nontrivial topology of the vortices, 
which carry an odd winding number of the order-parameter 
phase. These states—which are rather near-zero modes, 
because the topology may not impose a state at exactly zero 
energy in all cases—were predicted by the Bogoliubov-de 
Gennes equations^’^ and directly observed experimentally in 
NbSe 2 by scanning tunneling spectroscopy (STS), providing 
a striking verification of the theory. In cuprate high-T^ 
superconductors (HTS), similar STS experiments in vortex 
cores^~^^ failed to reveal the signature expected for bound 
states in a d-wave superconductor.This is surprising, 
because zero modes are in principle protected by topology 
and should be robust. Several explanations have been put 
forward,but a definitive interpretation of the vortex- 
core tunneling spectrum in high-T^, oxides remains to be 
found. By contrast, the iron-based high-T^ superconductors 
generally present vortex cores with the expected zero-energy 
peak characteristic of the bound states,^^"^^ although in one 
case there are core states but no peak.^^ 

In dirty superconductors, the zero-bias peak associated 
with the bound states is broadened,^^’^^ leading to a flat tun¬ 
neling spectrum in the core. This alone cannot resolve the 
cuprate puzzle, because the vortex-core spectrum is gapped 
at zero bias in these materials, with features reminiscent 
of the spectrum observed in the pseudogap phase above 
This has lead to the idea that the vortex cores are 
electronically similar to the mysterious pseudogap phase. It 
remains unclear how the pseudogap and superconductivity 


would interact in the vortex core and, in particular, how this 
interaction could release the topological frustration which 
demands a zero mode. 

Far from vortices, or in zero field, the low-temperature 
tunneling spectrum of bismuth-based HTS is rather well 
understood. It was recently found that an extension of the 
BCS theory, taking into account the band structure and a 
coupling to the antiferromagnetic spin resonance,^^ is able 
to reproduce the STS data of Bi 2 Sr 2 Ca 2 Cu 30 io +5 (Bi-2223) 
quantitatively.^^ This modeling shows that the interaction 
with the spin resonance changes the density of states sig¬ 
nificantly in zero field, by redistributing spectral weight 
over a broad energy range. A question naturally follows: 
how would this interaction change the electronic structure 
in a vortex? One possibility is that the antiferromagnetic 
order becomes static in the cores,^^ as several experiments 
have suggested.The local density of states (LDOS) in 
a vortex core with competing antiferromagnetic order may 
indeed share some similarities with the STS data for the 
cuprates.In the present work, we explore the oppo¬ 
site scenario, in which the antiferromagnetic fluctuations 
are not frozen, but remain dynamical in the vortex core. 
In contrast to the static case, the dynamical case cannot 
be formulated as a Hamiltonian mean-field problem: the 
coupling to the spin resonance enters via a nonlocal and 
energy-dependent self-energy. A simplified version of this 
model, ignoring nonlocal terms in the self-energy, was used 
earlier to study the effect of the spin resonance on the 
LDOS around a nonmagnetic impurity.Very recently, the 
same approach was applied in a small cluster to investigate 
charge-density wave formation.^^ Here, we solve this prob¬ 
lem for a vortex of d^ 2 _y 2 symmetry in a two-dimensional 
one-band system with parameters appropriate for Bi-2223. 
Although our results for the LDOS disagree with the vortex- 
core measurements in this material,^^ they show how the 
tunneling spectrum would look like in Bi-2223, in the ab- 
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sence of a pseudogap or on the strongly overdoped side, 
if the interpretation given in zero field in terms of Bogoli- 
ubov quasiparticles interacting with the spin resonance is 
correct. The spin fluctuations are a candidate for the pair¬ 
ing interaction in the iron-based superconductors, as they 
can generate the attraction for a pairing with symmetry 
between different Fermi-surface sheets in these multi-band 
systems. We believe that many of our results are relevant 
for the qualitative understanding of the recently measured 
vortex spectra in these materials. We will argue that the 
measurements reported in Ref. 31 contain signatures of the 
interaction with the spin resonance. 

The STS studies have so far been limited to—more pre¬ 
cisely, designed to—measuring locally, targeting the LDOS, 
which is the diagonal part of the real-space electron Green’s 
function. Some important properties of the quasiparticles, 
such as their nodal or antinodal character or their mean free 
path, do not appear clearly in the LDOS, but only indirectly, 
for instance via quasiparticle interference. After having 
discussed the LDOS, we will present another way of imag¬ 
ing the quasiparticles directly in real space. This approach 
probes the off-diagonal terms of the Green’s function and 
combines locality with nonlocal (momentum-space) sensi¬ 
tivity. Such a measurement is a considerable experimental 
challenge, requiring finely controlled double-tip tunneling, 
but could greatly enrich our knowledge of quasiparticles 
in correlated metals. We illustrate this by showing, in par¬ 
ticular, that the zero-energy quasiparticles remain nodal in 
the vortex core, despite the fact that the zero-energy LDOS 
extends along the antinodal directions. 

The model used and the methods employed to solve it are 
described in Sec. 11 and Appendices A and B. The results are 
presented in Sec. Ill: the self-consistent vortex order parame¬ 
ter in Sec. Ill A, the LDOS and a study of the spectral-weight 
transfer in Sec. IIIB, and the new quasiparticle imaging 
in Sec. Ill C and Appendix C. A discussion is proposed in 
Sec. IV. 


II. MODEL AND METHODS 

The Bogoliubov quasiparticles of a superconductor can 
emit or absorb spin fluctuations and thereby become short¬ 
lived, if they are coupled to the spectrum of spin excitations. 
In a translation-invariant superconductor, these inelastic 
processes are described by self-energy corrections in mo¬ 
mentum space, which modify both the normal (quasipar¬ 
ticle) and anomalous (gap) dispersions. At leading order, 
the self-energy is proportional to the convolution of the spin 
susceptibility with the quasiparticle propagator.^® In real 
space, the momentum convolution becomes a product and 
the self-energy takes the form 

£(r - r',icoJ = 

P !n„ 

x%csir ( 1 ) 


The self-energy is a matrix in Nambu space: the matrix ele¬ 
ment (II 22 ) describes the renormalization and lifetime 
of particlelike (holelike) quasiparticles and the matrix ele¬ 
ments II 12 and II 21 contain the renormalization and lifetime 
effects for the superconducting gap. The symmetry relations 
connecting these matrix elements are discussed below. In 
Eq. (1), the self-energy is a function of the fermionic Mat- 
subara frequency = (2n (5 with P = the 

inverse temperature, while the spin susceptibility Xs is a 
function of the bosonic frequency Q.^ = 2nn/P • ^BCS is the 
Green’s function describing the BCS-Bogoliubov quasiparti¬ 
cles in Nambu space, in the absence of coupling to the spin 
excitations. Finally, g is a coupling parameter. The justifi¬ 
cation for a constant (momentum-independent) coupling 
is that for the cuprates, in the energy range of interest, the 
spin susceptibility is dominated by the antiferromagnetic 
resonance,"^^ the weight of which is mostly localized near 
the momentum ( 71 , 71 ). 

Since the spin susceptibility has a sharp structure at the 
spin-resonance energy while the BCS Green’s function 
has structure at the gap edges ~ A, the main structure 
of the self-energy develops around + A, producing a 
kink in the quasiparticle dispersion and a dip in the tunnel¬ 
ing spectrum.®® The interplay of the van Hove singularity 
can induce additional structures and change these energies 
slightly. ®^’^® In the core of a BCS vortex, the main structure 
of the Green’s function comes from the bound states near 
zero energy. One may therefore expect that the scattering 
rate due to spin fluctuations is largest at the energy in 
the core, which would produce a dip at this energy in the 
final vortex-core spectrum. This naive expectation may miss 
part of the story, however, because it assumes a purely local 
effect of spin fluctuations, while in the homogeneous case 
the self-energy has a marked momentum dependence,®®’"^^ 
indicating significant nonlocal components. 

In order to compute the effect of spin fluctuations on 
the tunneling spectrum in a d-wave vortex, we replace the 
Green’s function %cs of a uniform BCS d-wave supercon¬ 
ductor in expression (1) by the BCS Green’s function 
calculated in the presence of a vortex. The latter is signifi¬ 
cantly modified with respect to ^BCS along with the forma¬ 
tion of the core states and, in particular, looses translation 
invariance: 

tir, r', ico„) = ^ g^Zs(r - r', in„) 

P in„ 

X r', iiOn - i^n)- (2) 

We neglect a possible feedback of the vortices on the spin 
susceptibility and assume that it remains translation invari¬ 
ant. We shall compute the self-energy (2) numerically on 
the real-frequency axis, as explained further below. An ex¬ 
ample of a numerical result is shown in Fig. 1. With the 
self-energy ready, the last step in order to obtain the vortex- 
core spectrum is to solve the modified Gorkov equations, 
written hereafter in matrix form with the Nambu indices 
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FIG. 1. Local self-energy at the center of an isolated vortex 
(solid red lines) and outside the vortex (dashed blue lines) for 
electronlike quasiparticles coupled to a (tt, n) mode of energy 
at T = 2 K. The dark gray regions indicates energies \oj\ 
the light gray regions correspond to < \cjo\ < + Aq. All model 

parameters are the same as in Sec. III. 


explicit and for an arbitrary complex energy 


-A-Si2(z) a 

f = (a 0^ 
(^21(2) ^ 22 (z)j (0 aj- 


(3) 


All products are matrix products in the implied spatial co¬ 
ordinates. = (^ + lJL)5rr' - trr' is 

the noninteracting normal-state Green’s function, w^ith jU 
the chemical potential and t^r' the hopping amplitude. 
A(r, rO is the superconducting pair potential describing a 
vortex with d^ 2_^2 symmetry. The symbol “T” means trans¬ 
position of the spatial coordinates and “t” means the same 
transposition followed by complex conjugation. It should 
be emphasized that the spin resonance is not the primary 
source of pairing in our model: the pair potential A(r, rO 
is generated by a different interaction of unspecified origin, 
via static Bogoliubov-de Gennes equations. We are primarily 
interested in the component ^ 11 , whose diagonal elements 
provide the local density of states: 

N(r,s) = - = siV). (4) 

The Dyson equation for is, as usual, ~ 

E, with the corresponding self-energy obtained by solving 
Eq. (3): 


S(z) = Sii(z)-[A + Si2(z)] 

X {[%-H-2)]^ + S22(z)}“'[A^ + S2i(z)]. (5) 

The difficulty of this problem stems from the combination 
of broken translational invariance and nonlocality in time. 
The rest of this section is mostly technical and describes our 
practical implementation of the solution. An overview of 
the calculation workflow is given at the end of the section. 

If the components E^^ of (2) are given, the calculation 
of the self-energy (5) for one particular energy requires a 
matrix inversion in the spatial indices. The quantity 


which is the renormalized propagator of the holelike quasi¬ 
particles, is indeed better written as (U — For 

practical reasons, such matrix inversions limit the system 
size to a few thousands lattice sites. One further matrix 
inversion is needed in order to obtain the local density of 
states, by solving the Dyson equation, also rewritten in 
the form ‘^11 = (U - %E) The reason for prefer¬ 

ring these rewritings, as opposed to simply solving, e.g., 
^11 = ^ — E) ^—which also requires only one ma¬ 

trix inversion because is known analytically—is that 
they allow for a better energy resolution by minimizing 
finite-size effects. Thanks to its translation invariance, 
the matrix % is best computed as a Fourier transform 
%ir,r',z) = - ?k) with the dis- 

persion measured from the chemical potential. We do this 
on a dense mesh of k vectors, with N = 1024 x 1024 points 
in the two-dimensional Brillouin zone, such that boundary 
effects are negligible. The expression (U — thus 

delivers a Green’s function free of boundary effects if E = 0, 
as opposed to the expression - E)“^. 

Let us now turn to the practical calculation of Eq. (2) 
for real frequencies. Following previous studies,^®’^^’"^® we 
use for the spin susceptibility a phenomenological model 
inspired by experiments: 


Xsir-r\inj = W,Fir-r') 


I — 

j -00 ^ 


( 6 ) 


1(e) represents a Lorentzian-broadened spin resonance at 
frequency and F(q) is peaked around the antiferromag¬ 
netic vector (n,n). In order to evaluate analytically the 
frequency sum in Eq. (2), we use the representation of the 
Matsubara Green’s function in terms of retarded (R) and 
advanced (A) functions 




_L r 

2n ico„ — s 

J -00 ^ 



p(.r,r',s) 
icOj^ — e 


(7) 


Performing the frequency sum and the analytic continuation 
^ CO + iO"^, we obtain the retarded self-energy on the 
real axis: 


t{r,r',co) = g^W,F{r -r') dsp{r,r',e) 




CO — e — e 


+ i0+ 


= a^F(r-rO dep(r,r',s)Bo(co,e), (8) 


(9) 


Boico,E) = -Bli-co,-E) 


= aM de[LrSe-^s)-LrSs + n,)] 

J —00 

/(-£) +b(s) 


( 10 ) 


CO — E — e + 10 + 
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The parameters o? and A are a dimensionless coupling 
and a typical energy scale (nearest-neighbor hopping), re¬ 
spectively (see Ref. 39). is the Lorentzian, with the 
energy width of the spin resonance, while f and b are the 
Fermi-Dirac and Bose-Einstein distribution functions, re¬ 
spectively. The function Bq(cl>,E) is analytically known and 
corresponds to the function B(co,E) of Ref. 39, in which 
the Dynes parameter F is set to 0"^. This function is peaked 
near E = co — El^sign(E), such that the energy integration in 
Eq. (8) is well convergent. The Dynes broadening is absent 
from the function Bq, because it is already implemented in 
the vortex Green’s functions 

We show now that the spectral functions Pij(r,r\e), 
needed for the calculation of the self-energy (8), can all 
be deduced from the element ‘11’ of the retarded Nambu 
matrix element ‘11’ denoted hereafter simply The 
vortex matrix Green’s function for a general complex energy 
2 is the solution of Eq. (3) with = 0. The element ‘11’ 
satisfies the Dyson equation with the 

BCS vortex self-energy given by setting = 0 in Eq. (5): 
Hvtx = For real frequencies, we have = 

3.nd 

- X! - ri, -CO - ir)A*(r', ra). (11) 

^1^2 

One can see (Appendix A) that the solution of Eq. (3) with 
= 0 has the analytic property ^uiz) = from 

which we deduce that G^^{e) = [G^^(e)]''' and further that 

Pii(r, r', e) = ^ s) - G*^[r', r, e)] (12) 

= p*i(r',r,e). 

Note that G^^(^r, r\ s) / Gytx(r^ r, e), such that the spectral 
function Pn is complex. We also have the property ^ 22(^3 — 
—‘^i*i(— 2 *), which implies 

P22(r, r', ^ [GvcxCr', r, -s) - G^^ir, r', -e)] (13) 

= Ph(r,r',-s). 

Together with the property (9), the relations (12) and (13) 
imply that the particlelike and holelike self-energies are 
related by Il 22 (r, r', co) = —E*^(r, r', —co). Finally, we 

have ^ 12(^3 = ^Ji(^*3 = ^21 (^“^*3 deduce 

P2i(r, r', ^ 04) 

= -P2iir', r,-s) = Pi2(r r,e), 
where, following the usual notation, the anomalous vortex 
Green’s function is F+^(£) = ‘^ 2 i(^ + i0^3- The latter can 
also be expressed in terms of Gy^x- Eq. (3) with = 0 gives 
(z), which means on the real axis: 

F+Jr,r',e) = 

- ^ %iri -r,-e- ir)A*(r 2 , r-^)G^{r 2 , r', e). (15) 


The relations (14) impose a connection between the pairing 
self-energies, namely r\ co) = E:*^(r, r\ -co). 

The numerical calculation runs as follows. 

1. Set the model parameters (dispersion pairing in¬ 
teraction) and determine the self-consistent BCS vor¬ 
tex pair potential A(r, rO; for this step, we use the 
Chebyshev expansion method described in Ref. 51 and 
briefly recalled in Appendix B. 

2. Choose the system size and use Eq. (11) and the cor¬ 
responding Dyson equation to calculate the vortex 
Green’s function Gy^x and to deduce with Eq. (15). 

3. Calculate and store the spectral functions pij using 
Eqs. (12-14) for all relevant energies and positions. 

4. Set the spin-resonance parameters, perform the en¬ 
ergy integration in Eq. (8) and store the self-energies 
So- 

5. Evaluate the modified vortex self-energy (5) on the 
real axis and solve the corresponding Dyson equation 
to obtain the Green’s function and finally the LDOS 
from Eq. (4) . 

All calculations of the self-energy (5) are performed at 
T = 2K, which is the base temperature of the experiments 
reported in Ref. 39. The self-consistent order parameter 
A(r, rO is computed at T = 0 for simplicity. 


III. RESULTS 

A. Self-consistent vortex pair potential 

A tight-binding model for the low-energy band struc¬ 
ture of Bi-2223 was obtained in Ref. 39, by fitting tun¬ 
neling spectra in zero field. We use the parameters 



r/a 


FIG. 2. Self-consistent vortex pair potential calculated for the 
nearest-neighbor bonds on the square lattice, using the Chebyshev 
expansion (lattice size: 1001 x 1001, expansion order: 2000; rela¬ 
tive convergence to self-consistency: < 6 x 10“^; see Appendix B). 
The color of dots in (a) and (b) varies from red (antinodal direc¬ 
tion) to blue (nodal direction). The dashed black and solid purple 
lines in (b) show two functional dependencies as indicated, with 
5 = 2.96a, <^o = 0.68a, and = 10.5a. 
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corresponding to the spectrum with a peak-to-peak gap 
of 44 meV These parameters are = 

(-206,56,-36,-10.3,27.9,-237) meV The bare BCS d- 
wave gap is Aq = 48.6 meV Figure 2 displays the self- 
consistent BCS pair potential A(r, r') calculated at T = 0 
for an isolated vortex. The pair potential vanishes in the core 
over a length scale similar to the bulk coherence length, con¬ 
sistently with previous Bogoliubov-de Gennes calculations 
for d-wave superconductors.The pair potential is 
nonzero only on the nearest-neighbor bonds (|r — r'| = a). 
The pairing strength was adjusted to reproduce the bulk 
gap Aq far from the vortex. The dots in Fig. 2(a) show 
the modulus of the pair potential for each bond; this rep¬ 
resentation differs from that in Ref. 52, where local d^ 2 _y 2 
and extended-s wave components were defined at each site. 
Beside the modulus shown in Fig. 2(a), the pair potential 
carries the signature, namely a plus (minus) sign 

on bonds running along x (y), as well as the topological 
phase given to an excellent approximation by —'d, where 
'd is the angle defined by the middle of the bond and the 
vortex center [see Fig. 2(a)]. 

The pair-potential modulus has cylindrical symmetry at 
large distances from the core, but presents some anisotropy 
at intermediate distances, as shown in Fig. 2(b). The gap 
increases faster along the (1,1) direction than along the 
(1,0) direction. By changing the band parameters, we have 
found that this behavior is model dependent. As already 
noted in Refs. 6 and 18, the relaxation of the gap is not very 
well described by the Ginzburg-Landau functional form, 
tanh(r/(^), whatever the value of The BCS expression of 
the coherence length is ^ = h{v^)/( tiAq). With our band pa¬ 
rameters, the average Fermi velocity is 2.63 x 10^ cm/s and 
thus ^ = 1.13 nm. The corresponding expected tanh(r/(^) 
dependency is shown in Fig. 2(b) as the dashed line. No 
good fit can be achieved with any value of We find that a 
different two-parameter functional form fits the numerical 
data much better, as indicated by the purple line. 

A self-consistent pair potential is nice to have, but not 
essential for the study of the LDOS in the core region.We 
have not found significant differences between the LDOS cal¬ 
culated using the Ginzburg-Landau and the self-consistent 
pair potentials. All data shown below correspond to the 
self-consistent case. 


B. Local density of states and spectral-weight transfer in 
the vortex coupled to the spin resonance 

The spin resonance is characterized by its energy = 
33.7 mey energy width = 4.2 mey and momentum width 
Aq = 1.15/a. With a coupling a = 0.7,^^ these parame¬ 
ters yield the zero-field DOS shown in the inset of Fig. 3, 
in very good agreement with the experimental tunneling 
spectrum.Note that the peak maxima in this zero-field 
DOS define a gap Ap = 46 mey slightly larger than the 
experimental gap of 44 meV We calculate the vortex-core 
LDOS in a 71 X 71 cluster having the vortex at its center.^^ 
In order to estimate the importance of finite-size effects, we 
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FIG. 3. Redistribution of spectral weight in the core of a d-wave 
vortex, due to interaction with an antiferromagnetic spin resonance 
of energy The black line with pink shade is the LDOS with the 
interaction turned on, while for the shaded blue line (BCS) it is 
turned off Vertical blue bars indicate the main spectral features of 
the BCS spectrum; pink bars indicate energies corresponding to the 
blue bars, shifted by ibH^. (Inset) Illustration of finite-size effects. 
Solid line: converged DOS in zero field. Dashed line: LDOS at the 
central site, after replacing the vortex by a uniform d-wave gap. 


replace the vortex pair potential by a uniform d-wave gap 
and compare the resulting LDOS at the central site with the 
fully converged DOS calculated in momentum space using 
1024 X 1024 k points (inset of Fig. 3). The good agreement 
between the two curves shows that the finite-size effects are 
small. Whenever this is possible, i.e., in the case without 
coupling to the spin mode, we also check our calculations 
against Chebyshev-expansion calculations in a much bigger 
cluster of size 1001 x 1001 (Appendix B). 

The LDOS curves calculated in the region of the core with 
and without the coupling to spin fluctuations are compared 
in Figs. 3-7. Figure 3 emphasizes the spectral differences 
at the vortex center. Fig. 4 compares two spectral traces 
along the directions (1,0) and (1,1), Figs. 5 and 6 illustrate 
the transfer of spectral weight, and Fig. 7 compares LDOS 
spatial maps at fixed energies. 

The interaction with the spin resonance leads to addi¬ 
tional structure in the spectra. The rearrangement of spec¬ 
tral weight can be qualitatively understood along the lines 
given in Ref 48: additional dips not present in the BCS 
spectrum correspond to energies where the scattering rate is 
large, due to enhanced emission or absorption of spin fluc¬ 
tuations. The additional peaks are less informative, because 
the spectral weight they carry is the one expelled from the 
dips towards both higher and lower energies. Schematically, 
the emission is strong at energies ^ if the BCS DOS at 
energy s — is large and, inversely, the absorption is strong 
ai £ < if the BCS DOS at ^ is large. The strong 
emission region marked as in Fig. 3 and characterized by 
a significant removal of spectral weight between ~ 20 and 
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60 mey as well as the absorption region marked aj, result 
from the zero-bias peak in the BCS spectrum: quasiparticles 
at these energies decay into the near-zero energy vortex-core 
states. Similarly the absorption and emission dips at a 2 and 
02 correspond to quasiparticles decaying into states at the 
gap edges, which survive as weak peaks at energies slightly 
lower than Aq in the BCS vortex-core spectrum. The peaks 
between and a 2 and between and 62 collect part of the 
spectral weight removed from the corresponding dips, but 
it is also seen that a great part of this weight is transferred 
to energies larger than ±100 meV The small peaks right at 
ai and seem more difficult to assess, but an inspection of 
the spectral traces in Fig. 4 reveals that they mark the onset 
of scattering at ±^ 5 . Below this threshold, the imaginary 
part of the self-energy (2) vanishes (see Fig. 1) and conse¬ 
quently the levels are renormalized to lower energy, but not 
broadened. In zero field or outside the vortex, the scattering 
rate tracks the linear increase of the d-wave BCS DOS and 
grows roughly linearly for |co| > the Kramers-Kronig 
related renormalization bends over with a weaker slope and 
as a result the coherence peaks develop a shoulder (white 
arrows in Fig. 4). In the vortex core, however, the scattering 
rate jumps at |co| = due to the zero-bias peak in the BCS 
spectrum and decreases for |co| > (Fig. 1); the renor¬ 
malization drops abruptly and changes sign, explaining the 
peaks at and e^. 

The spectral traces shown in Fig. 4 confirm the trends 
observed at the vortex center. The LDOS has more structure 
in the presence than in the absence of the coupling and 
dips in the interacting LDOS can be traced back to peaks in 
the BCS LDOS. On leaving the core, the dips at a 2 and 62 
become the usual dips of the zero-field spectrum (inset of 
Fig. 3). A tiny shift of these dips to lower binding energy in 
the core follows the tiny shift of the BCS coherence peaks. 
The coherence peaks themselves, which are washed out in 
the core due to both loss of superconducting coherence and 
increased scattering, start to develop near = 46 meV 
for r > 3a, where the modulus of the pair potential is 
close to its asymptotic value (Fig. 2). The small peaks at 
ai and become the shoulders on the coherence peaks as 
already discussed. Finally, the zero-bias peak, depleted form 
much of its spectral weight, splits into several structures 
reminiscent of the Caroli-de Gennes-Matricon bound states, 
which disperse with increasing distance from the center and 
whose intensity decreases over a length larger than the core 
size. Using the Chebyshev-expansion method on a large 
1001 X 1001 cluster, we have found that the wiggles in the 
BCS spectrum at energies above Aq in the core are finite- 
size artifacts: the converged spectrum is smooth at these 
energies. However, all little peaks at sub-gap energies in 
the BCS spectra of Fig. 4 are real. The energies and the 
amplitudes of these peaks depend on the band-structure 
parameters. A similar verification is not possible with the 
coupling turned on, but we believe that the sub-gap peaks 
in the interacting LDOS are real spectral features as well. 

It is interesting to study the spatial and energy depen¬ 
dence of the spectral-weight redistribution. To this end, we 
compare the spectral-weight transfer (SWT) induced by the 
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FIG. 4. Comparison of the LDOS in the vicinity of the vortex core 
with the coupling to the spin resonance turned off (left) and on 
(right). Upper and lower panels show a trace starting at the core 
center and going along the x axis and along the diagonal of the 
square lattice, respectively. The white and black arrows mark the 
energies ±^2^ and — (^2^ + Aq), respectively. 


spin resonance in zero field and around the vortex. In zero 
field, the SWT is defined in an energy domain as^^ 


1 

AW= — 
W 


dsN{s) 




:J— deJVBcs(e), 

"^BCS J\e\^s 


where W and are the total spectral weights, defined 
for our purposes as W = J|^|<; 2 oo variation 

of AW as a function of ^^lax ^ — [0, ^max] is displayed 
in Fig. 5(a) as the black line. Not much happens at sub-gap 
energies. The weight removed at the lowest energies is 
overcompensated, such that some weight is gained at 
and slightly above: this is the shoulder on the coherence 
peaks. The action really starts above A^, where 10% of the 
weight is removed over 50 meV and recovered at energies 
higher than 100 meV: this is the dip. In the vortex, we 
repeat the analysis by replacing the zero-field DOS by the 
average vortex LDOS over a square extending from — L to 
±L in both directions. There is a range of L values, between 
20a and 30a, where the result is almost independent of 
L [Fig. 5(a)]. We therefore set L = 25a in the following. 
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FIG. 5. Total (a) and local (b) spectral-weight transfer induced by 
the spin resonance (see text). The colored ranges in (a) correspond 
to the four energy domains in (b). Only half of the spatial region 
defined by L = 25a is shown in (b). The purple line represents the 
modulus of the pair potential as determined in Fig. 2. 

Smaller values of L focus too much on the core region, much 
larger values, if accessible, would mask the signature of the 
vortex and approach the zero-field result. With L = 25a, 
there is compensation of the average SWT very close to 
the two characteristic energies and A^. The question 
arises, whether some weight is redistributed spatially among 
neighboring sites: the answer is yes. In Fig. 5(b), the local 
SWT is plotted in four energy domains. The local SWT is 
defined as 

AW(r)=T dejv(r,e)-I dsNscs^r,e), 
J\e\^S l^BCS 

with W and Wbcs ^he average integrated spectral weights 
in the spatial region considered (L = 25a). At low energy 
1^1 < ^ 5 , a SWT of the order of 5% occurs from the vortex 
center to the surroundings. The depleted region is more 
extended than the core where the pair-potential modulus 
is significantly suppressed. The loss of low-energy weight 
in the vortex is obvious in Figs. 3 and 4, but the spatial 
compensation which occurs at distances r > 15a cannot 
be seen on those figures. Between and A^, there is 



(b) 



BCS BCS-l-spin resonance 


FIG. 6. (a) Total spectral-weight transfer induced by the vortex in 
the energy range \s\ < with and without coupling to the spin 
resonance. Three energy domains are defined in each case: below 
the maximum transfer (light-blue), between the maximum and the 
minimum (orange) and the high-energy region (red). The local 
spectral-weight transfers for the three domains are shown in (b). 

again a net SWT from inside the vortex to outside, although 
much weaker. The two small maxima at the vortex center 
in Fig. 5(b) are associated with the gain of weight near 
and Cl (see Fig. 3). Therefore, although there is no global 
SWT at the energy A^, locally there is a deficiency of weight 
in the vortex core and an excess outside. This unbalance 
is compensated in the region of the dip: between A^ and 
100 mey there is a loss of spectral weight everywhere, but 
much less in the core than outside, such that at 100 mey the 
global loss of weight is nearly uniform in space, as confirmed 
by the nearly uniform recovery above 100 meV 

From an experimental perspective, the SWT induced by 
the spin resonance is not accessible and thus not a quantity 
of particular interest: its measurement would require to 
turn off the interaction with the spin resonance. However, 
the SWT induced by the vortex itself is, in principle, directly 
accessible with nowadays technology, by switching the mag¬ 
netic field on and off and measuring the LDOS with and 
without the vortex. The local vortex-induced SWT is defined 
as 

AW(r)=T deJVo(e). 

J\e\eS ^ J\e\eS 

with ]V(r, e) the LDOS of the vortex and iVo(e) the zero-field 
DOS. The total vortex-induced SWT is obtained by replacing 
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FIG. 7. LDOS N(r,s) at several characteristic energies in the presence (bottom) and absence (top) of coupling to the spin resonance. 
The color scale goes from white (low) to dark blue (high) and covers the whole data range in each panel, such that absolute intensities in 
different panels cannot be compared.^^ 


N(r,s) by its average in the same spatial region as above 
(L = 25a). The result is displayed in Fig. 6 for the cases 
with and without the coupling to the spin resonance. There 
are qualitative differences, offering a chance to distinguish 
experimentally the two situations. The total SWT shows 
an accumulation of states at low energy [Fig. 6(a)]. Part 
of these states are outside the vortex core [Fig. 6(b)] and 
correspond to the increased DOS due to the Doppler shift of 
the dispersion;^® part of them correspond to the zero-bias 
peak in the core. Both the Doppler-shift and the core con¬ 
tributions are reduced by a factor close to two when the 
interaction in present. This is consistent with a renormaliza¬ 
tion of the quasiparticle velocity by a factor 1/(1 + A), since 
the renormalization factor A is close to unity.®^ The Doppler- 
shift approximation is valid at low energy and indeed one 
sees that it breaks down when approaching the gap scale. 
In this energy domain, the LDOS outside the vortex is lower 
than in the absence of field and as a result the total SWT 
decreases. In the BCS case, the compensation is complete 
at the energy Aq. This is only true on average, though, 
because locally there is still excess weight in the core and 
deficiency outside at the energy Aq. This is restored by a 
very nonuniform SWT at higher energy, as seen in Fig. 6(b). 
Here we observe the most significant differences. With the 
coupling on, there is no compensation on average at the 
peak energy A^ and the local compensation occurring near 
60 meV at the minimum of the total SWT is almost com¬ 
plete, such that the high-energy SWT is small and much 
more uniform spatially than in the BCS case. 

We close this section by comparing in Fig. 7 the spatial 
LDOS distributions in the presence and in the absence of the 
coupling, at a few characteristic positive energies. Similar 
behaviors are found at negative energies. The most signif¬ 
icant differences occur at ^ and e = + Aq. The 

zero-energy states are similar in both cases: they are local¬ 
ized on a few lattice sites and spread along the principal 
axes, not along the nodal directions. At ^ the BCS 


case shows a square, which resembles a Caroli-de Gennes- 
Matricon state. The radius of such a state at energy E may 
be estimated as®^ = (E/Ao)[l + 2r^/(n^)]r^, where 
represents the core radius. For E = this expression 
gives the value ~ 5a indicated by the data if one takes 
= 3.9a, which is a reasonable value, considering the gap 
profile shown in Fig. 2. The spatial structure of the LDOS 
is completely changed by the spin resonance: the density is 
spread out of the vortex (see also Fig. 5), such that there 
is no (quasi-) localized state at this energy. At the gap edge, 
the LDOS shows in both cases a depression in the core, but 
the latter appears to be shrunk and rotated by 45° by the 
spin resonance. The difference is most spectacular at the 
dip energy + Aq. In the BCS case, there is a deficiency of 
density in the core (see also Fig. 6) and the LDOS displays 
quasiparticle interference patterns. The latter are due to 
scattering on the vortex, but also on the cluster’s bound¬ 
aries. We have checked this by repeating the calculation in a 
much bigger cluster with the Chebyshev-expansion method. 
These interference patterns disperse with energy, as seen 
in the 150 meV map. They are also present a lower en¬ 
ergies, but too weak to be clearly resolved in comparison 
with the vortex-related structures. With the spin resonance, 
there appears to be an excess rather than a deficiency of 
density in the core. As already discussed, + Aq is the 
energy where the scattering rate is largest outside the vor¬ 
tex, due to the BCS coherence peaks at Aq. Although also 
large, the scattering rate is about two times smaller in the 
core where these coherence peaks are reduced (see Fig. 1). 
Hence less density is removed in the core than outside. It 
should be noted that this difference is small, of the order 
of 6%, but is magnified by the choice of the color scale 
in Fig. 7. For comparison, the LDOS difference at ^ = 0 
between the maximum and the minimum is 90%.®^ Even 
more striking at + Aq is the complete washing out of the 
interference patterns by the spin resonance. This suggests 
that the mean free path does not exceed the quasiparticle 
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FIG. 8. Nonlocal conductance ratio at zero energy in zero field (columns 1 and 2) and close to a vortex (columns 3 and 4). is fixed 
to (0,0), corresponding to the center of the panel, in columns 1, 2, and 4; = (— 20,5)a in column 3 (black dot). The yellow circles 

indicate the vortex center. All panels in a row share the same color scale displayed on the right. 


wavelength. We may tentatively estimate the mean free 
path as £ = (v*)t, where (v*) is the average renormalized 
Fermi velocity, which is 1.15 x 10^ cm/s, t = /i/(2Zy) is 
the quasiparticle lifetime, y being the scattering rate, and 
^ the quasiparticle residue. Outside the vortex, 

y is typically 100 meV at the energy + Aq and 50 meV 
at higher energies (Fig. 1). The corresponding values of £ 
are 2.3a and 4.5a, respectively. The former is shorter than 
the period of the BCS quasiparticle oscillations in Fig. 7, 
while the latter is comparable. This may explain why no 
interference at all is seen at + Aq, while some traces 
remain at 150 meV Our naive formula underestimates the 
mean free path, however: a more rigorous evaluation, to be 
performed in the next section, gives values nearly five times 
larger. 


C. Bogoliubov quasiparticles in real space 

The quasiparticle LDOS discussed in the previous section 
can be measured using a scanning tunneling microscope 
(STM).^^ There are other important quasiparticle proper¬ 
ties which are not visible in the LDOS alone. For instance, 
their nodal character: the fact that the zero-energy LDOS in 
Fig. 7 extends along the lattice axes rather than along the 
diagonal does not imply that the zero-energy quasiparticles 
in the vortex have lost their nodal character. This nodal 
character is still present, as we will see. As another example, 
the short mean free path of the quasiparticles coupled to 
the spin resonance is not directly apparent in the LDOS, but 
only indirectly through its effect on the interference pat¬ 
terns. It is therefore useful to go beyond the LDOS and to 
consider nonlocal effects that reveal directly these quasipar¬ 
ticle properties. While the LDOS is encoded in the diagonal 
elements G(r,r,s) of the retarded Green’s function G, the 


off-diagonal elements G(r,r^, ^) contain all additional in¬ 
formation about the single-particle excitations. It is well 
established that these off-diagonal terms are accessible ex¬ 
perimentally, for instance by local double-tip tunneling. 

A double-tip STM with nanometer distances between the 
tips has already been demonstrated.^^ Here, we consider 
the quantity 


Im[G(ri, r 2 , g) + G(r 2 , r^, g)] 
Im[G(ri, ri, s) + G(r 2 , r 2 , s)] ’ 


(16) 


We show in Appendix C that this quantity is the relative 
change of the tunneling conductance measured at zero tem¬ 
perature with two coupled tips at points and r 2 , with 
respect to the conductance measured with two uncoupled 
tips at the same two positions. We refer to this as the non¬ 
local conductance ratio (NLCR). In the proposed setup, the 
two tips are coherently connected to the same reservoir. As 
a result, the NLCR can be positive or negative, unlike in the 
setup of Ref. 60, where the “transconductance” is positive 
definite. In a translation-invariant system, R(ri — r 2 , s) is 
proportional to the imaginary part of the Fourier transform 
of the Green’s function G(k, s). Because G(k, e) = G(-k, e), 
this is identical to the Fourier transform of the imaginary 
part of G(k, e), i.e., the Fourier transform of the spectral 
function A(k,E). For ^ < 0, this quantity could in principle 
be deduced from photoemission data, but this has not been 
reported so far, to the best of our knowledge. 

The NLCR at zero energy, calculated without magnetic 
field as well as close to a vortex, are compared in Fig. 8. 
Since R(ri, r 2 , s) has a geometrical 1/r decay for long dis¬ 
tances, we multiply by | — r 2 1 for plotting. In zero field, 

the NLCR is shown for large distances (first panel in each 
row) and for short distances over the area corresponding 
to the cluster size of the vortex calculation (second panel 
in each row). The nodal character and the long-range na¬ 
ture of the zero-energy excitations appear very clearly in 
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Zero field 


R(ri, r2, e = n, + Ao) X |ri - r2l 

Vortex core 



FIG. 9. Same as Fig. 8, at the energy e = + Aq. The dashed circle indicates the mean free path, calculated by fitting an exponential 

decay to the NLCR. 


this representation, both without and with the coupling to 
the spin resonance. The period of oscillation corresponds 
to the nodal Fermi wavelength, as indicated in the figure. 
In the BCS case, the Fermi wavelength is almost exactly 
commensurate with the lattice, 2n/k^^ = 3.003(a,a). It is 
renormalized to slightly lower values 2n/k*^ = 2.961(a, a) 
by the spin resonance: this small incommensurability ex¬ 
plains the Moire pattern with a period ~ 36a along the 
diagonals. 

In order to visualize excitations in the vicinity of the vor¬ 
tex core, we first take = (-20,5)a and plot R(_ri, r 2 ,0) 
as a function of r 2 . This may be loosely understood as the 
long-time behavior of the wave function for a particle cre¬ 
ated at r^. It is seen that the latter is similar to the zero-field 
quasiparticle. The spreading away from the nodal directions 
appears to be wider than in the zero-field case, but this is 
a finite-size effect. The scattering on the vortex is clearly 
visible; one also sees that it is a relatively weak effect. Note 
that the Doppler effect, by which all momenta are shifted 
by a quantity (h/2r)e^, where r is the distance to the vortex 
and is the direction of the supercurrent, is too small to 
be seen. For r ~ 20a, the relative change of wavelength is 
1/(1 + 2A:pr) ~ 2%. When the particle is created right at the 
vortex center (rightmost panel in each row), the NLCR is 
globally smaller because the LDOS at this point is large. The 
excitation remains nevertheless nodal, despite some loss of 
nodal intensity (which is less pronounced in the presence of 
the spin resonance), in contrast to the LDOS, which extends 
along the antinodal directions at ^ = 0 (Fig. 7). 

Figure 9 shows the NLCR calculated at the energy ^ H- 
Aq, where the effects of the spin resonance are strongest. In 
the BCS case, the spectral weight is mostly in the antinodal 
region, such that the NLCR extends mainly along (1,0) and 
(0,1). At this energy, the antinodal wave-vector is = 
(0.22, l)7r/a, corresponding to a period of 9.1a along one 
direction and 2a along the other, as indicated in the figure. 
Some nodal states are also mixed in, but their contribution 


is weak. Note that the scattering on the vortex is inexistent 
at this energy, both in the BCS and spin-resonance cases. 
Three major differences are seen with the coupling to the 
spin resonance. First, the amplitude of nonlocal effects is 
reduced by an order of magnitude (see color scales on the 
right of Fig. 9). Second, the (tt, tt) scattering suppresses 
spectral weight and broadens states in the antinodal region, 
such that the better-defined nodal states are dominating the 
shape of the NLCR. Third, the appearance of a mean free 
path is manifested by an exponential decay of the NLCR 
with increasing distance Ir^ — r 2 |. Fitting this exponential 
decay, we obtain a mean free path I = 11a (dashed circle in 
the figure), which is much larger than our previous estimate 
based on the average Fermi velocity and scattering rate. 
Repeating the same fitting as a function of energy, we find 
that the mean free path tracks the energy dependence of the 
average scattering rate, but with an exponent different from 
— 1. Both quantities are approximately related by £/a 
320 X (y/meV)“^/^. In this analysis, the average scattering 
rate was defined as ^(e) = (l/2)(p(k,E)Im[-Eii(k,E) - 
^)])bz. with p(k, e) = Ao(k, ^)/(Ao(k, s))bz, Ao(k, e) 
being the spectral function without pairing gap (A^ = 0) 
and (• • • )bz standing for a Brillouin-zone average. 


IV DISCUSSION 

The model (1) reproduces the zero-field STM data of 
Bi-2223 very well;^^ yet, the same model (2) fails to re¬ 
produce the experimental vortex-core spectrum of this 
material.^^ The latter is similar to published results for Bi- 
2212,^^"^^’^^’^^ lacking a zero-bias peak in the core, but 
showing weak low-energy structures on top of a (pseudo- 
) gapped spectrum. The self-energy (2) does not suppress 
the zero-bias peak, although it reduces its weight; neither 
does it smear out the BCS spectra in a way which would 
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turn the vortex-core spectrum into the structureless signa¬ 
ture observed in dirty superconductors, but rather it adds 
structure to it. Lastly while the experiment suggests that 
the dip is disappearing in the vortex core, it is not the case 
with the model (2). The fact that the model (2) yields a 
peak in the vortex core illustrates the robustness of this zero- 
bias feature, associated with the topological defect carried 
by the vortex. It is unlikely that any theory in which the 
zero-field spectrum is entirely made of (possibly damped) 
Bogoliubov quasiparticles of a BCS superconductor can pro¬ 
duce a vortex-core LDOS without zero-bias peak, unless 
something that goes beyond the BCS theory happens in 
the vortices. The nucleation of a static antiferromagnetic 
order is one possibility,^^’^^’^^ but a detailed comparison 
of this theory with STM data of the cuprates is still miss¬ 
ing. Charge order in the vortices has also been suggested 
experimentally,^^ but the local spectral signatures of such 
an order remain so far unknown. While these interpreta¬ 
tions focus on the possible non-BCS nature of the vortices, 
many experiments point to a non-BCS nature of the zero- 
field spectrum as well, in relation to the phenomenon of the 
pseudogap. In this perspective, the success of the model 
(1) in reproducing the zero-field data must be an accident, 
indicating that the general BCS-Eliashberg structure of the 
Green’s function would be correct, despite the fact that the 
underlying physics would be wrong. The failure of (2) in the 
vortex core would then reveal the trick, because there the 
BCS-Eliashberg structure, which assumes superconducting 
coherence of all excitations, would be inappropriate. Dis¬ 
criminating between the two scenario requires a complete 
microscopic theory of the pseudogap and its interaction with 
superconductivity, which is not yet available. 

In ferropnictides, the role of antiferromagnetic spin fluc¬ 
tuations as the driving force for pairing is more firmly estab¬ 
lished than in the cuprates.Tunneling spectra in zero field 
and in vortices are available for a number of these materi¬ 
als and spectral structures suggestive of a coupling to the 
spin resonance, as measured by neutron scattering, are rou¬ 
tinely observed.For example in Bao. 6 Ko. 4 ^^^ 2 As 2 , 
the (tt, n) resonance is measured at 14 meV,^^ the tunneling 
spectrum shows a gap of 7 meV and presents the expected 
DOS structure near 21 meV.^^ The evidence of the coupling 
in the vortex-core spectra is more elusive, but encouraging. 
In FeSe, for instance, the vortex-core spectrum reported in 
Fig. 2A of Ref. 31 has a lot in common with that displayed 
in Fig. 3. The weak structures at ±3 meV in the former 
are reminiscent of the ones observed at the onset of scat¬ 
tering in the latter (a^ and e^); the structures at ±5 meV 
correspond to the main absorption and emission at a 2 and 
02 . These numbers are consistent with the value of the 
zero-field gap, observed near 2 meV, and the value of the 
spin resonance energy which, although not yet measured 
by neutron scattering, is expected near 4Ak^Tc = 3 meV 
in this material.Further evidence could emerge from a 
careful study of the spatial structure of the LDOS, or from 
an analysis of the vortex-induced spectral weight transfer. 

The study of quasiparticle interference (QI) by STM has 
been a fruitful development in the last decade. QI does 


not image the quasiparticles directly, but their interferences 
due to multiple scattering on defects: this is an advantage, 
because the interferences are long-ranged even if the quasi¬ 
particles are not, but it is also a weakness of the method, 
which makes it impractical for clean systems. Measuring the 
nonlocal conductance ratio (NLCR) introduced in Sec. Ill C 
would provide additional information about the quasiparti¬ 
cles. This is a considerable challenge by tunneling, perhaps 
less so by photoemission, although with the latter technique 
the local information is lost and positive energies are not ac¬ 
cessible. The NLCR exists also in clean systems and images 
the quasiparticles directly, showing their momenta—^while 
QI shows differences of momenta—as well as their spatial 
coherence range. Via the NLCR, local tunneling can probe 
the reciprocal space with two advantages compared to pho¬ 
toemission: an easy access to positive energies and the 
possibility to follow local variations of the reciprocal-space 
properties. 

V CONCLUSIONS 

We have calculated the electronic properties of a vortex 
in a two-dimensional one-band d-wave superconductor, in 
which the Bogoliubov quasiparticles interact with a spin 
resonance centered at momentum ( 71 , 71 ). Unlike previous 
calculations dealing with a static antiferromagnetic order 
in the vortex core, the dynamical case considered here is 
not amenable to a mean-field treatment, but requires to 
evaluate a self-energy that is nonlocal in space and time. 
We have discussed several signatures of the coupling to the 
spin resonance in the local density of states near the core. 
In spite of the fact that the model we use fits quantitatively 
the tunneling spectrum of Bi 2 Sr 2 Ca 2 Cu 30 io +5 in zero field, 
it fails to reproduce the peculiar vortex-core spectra of the 
cuprates. We believe that passing the test of the vortex- 
core spectrum is a tough sanity check for all theories of the 
cuprates electronic structure. Our results may nevertheless 
be useful to understand the vortex-core tunneling spectra of 
iron-based superconductors and we have argued that some 
of the signatures we discuss may be present in published 
data for FeSe. 

Bogoliubov quasiparticles coupled to spin fluctuations 
loose spatial coherence. We have shown that this is mani¬ 
fested by an extinction of quasiparticle interference at the 
energies where the effect of the coupling is strongest. We 
have also discussed a new way to look at quasiparticles in 
real space, which allows one to access quasiparticle proper¬ 
ties that are not directly visible in the LDOS, such as their 
nodal/antinodal character and their mean free path. 
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Appendix A: Symmetry properties of the Green’s function 

In the absence of coupling to spin fluctuations = 0), 
the Nambu-Gorkov equation (3) reads 

-A ) f ^n(^) _ fn 0) 

[ -A^ t^2i(z) ^22(z)J [o Ti-j’ 

with the explicit solution 

^21 (z) = -{[%’' (-z)A%o(z)]"' + A}-' 

^ 22 (z) = + A^%{z)A}-\ 

The property %(_z) = easily checked from the 

Fourier representation of %, implies that 

= ^11 (^*), ^ 12 (^) = ^ 22 (^) = ^ 22 (^*)* 

Furthermore, the property A = A^, which follows from the 
S)mimetry of the pairing interaction, allows one to deduce 
two additional relations: 


^22(^) = -^n(-^*), ^12W=^2 i(-^*)- 


Appendix B: Chebyshev expansion of the Green’s functions 


Consider a one-band quadratic Hamiltonian defined on 
a lattice spanned by the discrete vectors r. In the super¬ 
conducting state, there are two degrees of freedom at each 
lattice site, namely the Bogoliubov-de Gennes amplitudes 
[u(r),v(r)]. The real-space retarded Green’s function at 
energy E is G(r, r\E) = (r |[(E + iO"^)Il In this 

expression, H is the matrix representing the Hamiltonian 
and the notation | r) means a state describing an electron 
localized at point r: the corresponding state vector has 
u(r) = 1 and all other components equal to zero. For each 
pair (r, rO, the Hamiltonian is the 2x2 block 


'-/xSrr' + trr' A(r,rO 
^ A*(r» 


(Bl) 


The expansion method of Ref. 51 consists in expanding 
the matrix G(E) = [(E + i0"^)ll - on Chebyshev poly¬ 
nomials in H. The Chebyshev polynomials are r„(x) = 
cos(narccosx), defined in the range -1 ^ x ^ 1. The 
expansion of G(E) on these polynomials reads 


G(E)=-V ,- - 


n=0 


Vl -£2 


(B2) 


Tildes indicate rescaled dimensionless energies, which fall 
within the range [—1,1] on which the Chebyshev polyno¬ 
mials are defined: H = (H — b)/a, E = (E — b)/a, with a 


and b the width and center of the spectrum of H, respec¬ 
tively (exact values are not required). The polynomials obey 
the recursion relation T^+iCx) = 2xr„(x) — r„_i(x), such 
that the evaluation of the matrix elements (r |r„(H)|rO can 
be performed iteratively. This computational method is 
very appealing for several reasons. First, the calculation 
of the matrix elements only requires to evaluate matrix- 
vector products Hlip), with no need to store the Hamil¬ 
tonian in memory; only three state vectors need to be 
stored in the recursive scheme: — kO. IV^i) = 

IV^n+i) — — |'0n-i)- Second, when the matrix ele¬ 

ments are known, the Green’s function can be computed at 
any energy E in almost no time. Third, the calculation is triv¬ 
ially parallel in the positions r. Lastly, as the Hamiltonian 
propagates the initial state \r') on the neighboring sites of 
r' and so on at each iteration, the linear lattice size needed 
to calculate the matrix element {r\H^\r') is proportional 
to n. Nevertheless, manageable lattice sizes give accurate 
results thanks to the good convergence properties of the 
Chebyshev expansion. 

When the n sum in Eq. (B2) is truncated to some maxi¬ 
mal value N, the linear lattice size M must ideally be such 
that boundary effects do not affect the last matrix element 
{r\Tj^(H)\r^). For instance, if the propagation proceeds only 
through hopping to the nearest neighbors, it takes N = 2M 
iterations until the reflection from the boundary propagates 
back to the site r and N = 4M iterations until interferences 
between the reflections on the two opposite boundaries can 
be felt at r (assuming that r is the central site of the system, 
which can always be arranged). If the minimal size re¬ 
quirements are met, there remain nevertheless oscillations 
due to the truncation itself, known as Gibbs oscillations. 
Those can be suppressed by well-known procedures.^® We 
use the Lorentz kernel, which amounts to multiplying each 
term of the sum in Eq. (B2) by sinh(fiV — fn)/sinh(fiV) 
and is equivalent to introducing a phenomenological Dynes 
scattering rate F = af. 

The BCS gap equation is A(r,rO = —y(r, rO(Cr|Cr/j), 
where V(r,r') is the pairing interaction and annihi¬ 
lates a spin-cj electron at position r. The average value 
{Cr^Cr'i) can be related to the retarded anomalous function 
F^(r,r\E) = {r\G(E)\r'). \r) is the state vector describing 
one hole localized at point r, i.e., with v(r) = 1 and all 
other components equal to zero. We obtain 


A(r,rO = nr,rO dE/(E) 


X [£+(r,r',£)-£+(r',r,-£)]*, (B3) 

Zn 


where /(E) = + i) i is the Fermi function. Setting 

the temperature to zero, inserting the Chebyshev expansion 
of E"^ into Eq. (B3) and performing analytically the integral. 
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we are led to the following expression: 


A{r,r') = V{r,r')-y] 


_oo _ ( g-marccos(-b/a) 


X [(r|T„(H)|r') + (f'|T„(H)|r)] I . (B4) 


The term of order n = 0 disappears because (r|ro(H)|rO = 
(r |r0 = 0. Equation (B4) must be solved self-consistently, 
with H given by Eq. (Bl). 


Appendix C: Double-tip tunneling 

Consider the usual electron tunneling problem, with two 
electrodes spanned by the vectors I (left) and r (right) 
and characterized by the single-particle retarded Green’s 
functions G(l,V,s) and G(r,r\s), respectively The two 
electrodes are coupled by a tunneling matrix element t(Z, r) 
and a chemical potential difference eV is applied. The 
single-particle current at leading order in t(Z, r) is^^ 

2716 r °° 

ds[f{s-eV)-f{s)] 

J —00 

X (Cl) 

IVrr' 


p(r,r',s) is the spectral function, which is related to the 
retarded Green’s functions in the same way as in Eq. (12). 
Although this formula was initially derived for quadratic 
Hamiltonians in the electrodes, it can be shown that the 
insertion of the interacting spectral functions in Eq. (Cl) 
correctly accounts for all correlations present in the elec¬ 
trodes, only neglecting correlations that are induced by the 
tunneling term.^^ In the proposed setup (Fig. 10), the tun¬ 
neling matrix element vanishes unless I = Iq and r = or 


r = r 2 . We can write We 

insert this into Eq. (Cl) and note that p{Iq, Iq, s -eV)is the 
LDOS in the left electrode representing the tip. We adopt the 
standard approximation and take it as a constant. Likewise, 
p(r,r,a) = (-l/7r)ImG(r,r,a) = iV(r,a) is the LDOS in 
the right electrode. Since the bias V only appears in one 



FIG. 10. Principle of double-tip tunneling to a common reservoir. 
The whole double-tip system (orange) is assumed to have the same 
density of states. 

of the Fermi functions, we can differentiate and obtain the 
tunneling conductance as 

di r °° 

— oc da[-/'(s-eV)] {t^JV(ri,e) + t^JV(r 2 ,e) 

J -00 

+tit2 (-^)lm[G(ri,r2,s) + G(r2,ri,e)]}. (C2) 

The first two terms in the curly braces give the parallel 
conductances associated with the two tips, proportional to 
the thermally-broadened LDOS at the points and r 2 . The 
third term accounts for transport processes involving both 
tips. The latter is sensitive to the quasiparticle phase change 
between the points and r 2 and thus provides a correction 
to the conductance which can be positive or negative. If the 
amplitudes and t 2 are equal, it is seen that the relative 
change of conductance induced by the two-tip processes is 
given at zero temperature and for the energy s = eV by 
Eq. (16). 
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